МІНIСТЕРСТВО ОСВIТИ І НАУКИ УКРАЇНИ
НАЦІОНАЛЬНИЙ УНІВЕРСИТЕТ “ЛЬВІВСЬКА ПОЛІТЕХНІКА”
РОЗРАХУНКОВО-ГРАФІЧНА РОБОТА №2
ЗНАХОДЖЕННЯ ДИНАМІЧНИХ ХАРАКТЕРИСТИК
ТЕХНОЛОГІЧНОГО ОБ’ЄКТА ШЛЯХОМ РОЗВ’ЯЗУВАННЯ
СИСТЕМИ ДИФЕРЕНЦІАЛЬНИХ РІВНЯНЬ
Львів - 2008
Для застосування схеми числового інтегрування до системи звичайних диференціальних рівнянь (ЗДР) необхідно її розв'язати відносно перших похідних шуканих функцій. В літературі з числових методів традиційно формули числового інтегрування системи двох ЗДР методом Рунге-Кутта записують у вигляді:
де , - праві частини ЗДР, h - рівномірний крок вибору вузлів інтегрування.
Завдання
Для протічної гідравлічної ємності з коротким та довгим трубопроводами відповідно на вході та виході (рис.1) математична модель (1) має вигляд:
(1)
де , , , , .
Знайти розв’язок системи звичайних диференціальних рівнянь (ЗДР) - H(t), Q1(t) (зміну рівня в ємності та витрати рідини в трубопроводі) на інтервалі t([0; 4], c з кроком Δt=0,05 c, якщо
d=0,25 м
r1=0,06 м
L1=90 м
r2=0,06 м
L2=2 м
P1=1500 Па
P2=300 Па
(=0,9
ρ=1e3 кг/м3
g=9,81 м/с2
ν=1e-5 Па
Початкові умови: значення рівня H(0)=0,1 м та витрати Q1(0)=0,001 м3/с.
Розв’язок подати у вигляді графіка.
Розв’язання
Застосуємо цю схему числового інтегрування системи двох ЗДР для системи (1) диференціальних рівнянь. Спочатку розв’яжемо її відносно перших похідних шуканих функцій, тобто
(2)
Позначимо
;
або
.
Тоді схема числового інтегрування системи двох ЗДР (2) матиме вигляд: (3)
Програма числового розв'язування системи двох ЗДР мовою С.
#include <stdio.h>
#include <conio.h>
#include <math.h>
float f1(float t, float H, float Q1)
{float d=0.25,r2=0.06,L2=2,g=9.81,ro=1e3,nu=1e-5;
float S=M_PI*d*d/4,P2=300;
float k2=M_PI*pow(r2,4)/8*nu*L2;
return 1/S*(Q1-(k2*ro*g*H-P2)/ro);
}
float f2(float t, float H, float Q1)
{float r1=0.06,L1=90,dz=0.9,ro=1e3;
float k1=2*M_PI*r1*r1*sqrt(r1/L1*dz), P1=1500;
float A=4*M_PI*pow(r1,3)/dz;
return 1/A*(k1*k1*P1/ro-Q1*Q1);
}
main()
{clrscr();
float t, t0=0, tf=4, n=40;
float h,k1,k2,k3,k4,l1,l2,l3,l4,Hi=0.1,Q1i=0.001;
h=(tf-t0)/n;
for(t=t0;t<=tf;t+=h)
{printf("\n t=%f Hi=%f Q1i=%f",t,Hi,Q1i);
k1=h*f1(t, Hi,Q1i);
l1=h*f2(t, Hi,Q1i);
k2=h*f1(t+h/2, Hi+k1/2, Q1i+l1/2);
l2=h*f2(t+h/2, Hi+k1/2, Q1i+l1/2);
k3=h*f1(t+h/2, Hi+k2/2, Q1i+l2/2);
l3=h*f2(t+h/2, Hi+k2/2, Q1i+l2/2);
k4=h*f1(t+h, Hi+k3, Q1i+l3);
l4=h*f2(t+h, Hi+k3, Q1i+l3);
Hi=Hi+1.0/6*(k1+2*k2+2*k3+k4);
Q1i=Q1i+1.0/6*(l1+2*l2+2*l3+l4);
}
getchar();
return 0;
}
Запишемо в таблицю 1 кожне десяте значення функцій для рівня H(t), м та витрати Q1(t), м3/с.
Таблиця 1
t,c
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
H(t),м
0.1000
0.7132
1.3263
1.9394
2.5525
3.1656
3.7786
4.3916
5.0045
5.6175
Q1(t), м3/с
0.0010
0.0010
0.0010
0.0009
0.0009
0.0009
0.0009
0.0009
0.0009
0.0009
t,c
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
H(t),м
6.2304
6.8433
7.4562
8.0691
8.6820
9.2948
9.9076
10.5205
11.1333
11.7461
Q1(t), м3/с
0.0009
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
t,c
2
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
H(t),м
12.3588
12.9716
13.5843
14.1971
14.8098
15.4226
16.0353
16.6480
17.2607
17.8734
Q1(t), м3/с
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0007
t,c
3
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
4
H(t),м
18.4861
19.0987
19.7114
20.3241
20.9367
21.5494
22.1620
22.7747
23.3873
23.9999
24.6126
Q1(t), м3/с
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
Розв’яжемо систему двох ЗДР засобами MatLab.
Створимо файл dataf.m, в якому запишемо вхідні дані та обчислимо коефіцієнти.
%data5
ro=1e3; mu=1e-5; g=9.81; dz=0.9;
d=0.25;
S=pi*d^2...
